arXiv:1504.07907v2 [cs.CV] 30 Apr 2015 


A Flexible Tensor Block Coordinate Ascent Scheme for Hypergraph Matching 


Quynh Nguyen Ngoc^’^, Antoine Gautier^, Matthias Hein^ 

^Max Planck Institute for Informatics, Saarbriicken, Germany 
^Saarland University, Saarbriicken, Germany 


Abstract 

The estimation of correspondences between two images 
resp. point sets is a core problem in computer vision. One 
way to formulate the problem is graph matching leading to 
the quadratic assignment problem which is NP-hard. Sev¬ 
eral so called second order methods have been proposed 
to solve this problem. In recent years hypergraph match¬ 
ing leading to a third order problem became popular as it 
allows for better integration of geometric information. For 
most of these third order algorithms no theoretical guaran¬ 
tees are known. In this paper we propose a general frame¬ 
work for tensor block coordinate ascent methods for hyper¬ 
graph matching. We propose two algorithms which both 
come along with the guarantee of monotonic ascent in the 
matching score on the set of discrete assignment matrices. 
In the experiments we show that our new algorithms outper¬ 
form previous work both in terms of achieving better match¬ 
ing scores and matching accuracy. This holds in particular 
for very challenging settings where one has a high number 
of outliers and other forms of noise. 


1. Introduction 

Graph resp. hypergraph matching has been used in a va¬ 
riety of problems in computer vision such as object recogni¬ 
tion [20], feature correspondences [7, 24], shape matching 
[1 1, 23, 25] and surface registration [27]. Given two sets of 
points, the task is to find the correspondences between them 
based on extracted features and/or geometric information. 
In general, the graph matching problem is NP-hard, there¬ 
fore, many approximation algorithms have been proposed 
over the years. They can be grouped into second order and 
third order methods. 

Among recent second order approaches, Leordeanu and 
Hebert [17] proposed the Spectral Matching (SM) algo¬ 
rithm, and Cour et al. [9] the Spectral Matching with Affine 
Constraint (SMAC). Both of these methods are based on the 


best rank-1 approximation of an affinity matrix. Leordeanu 
et al. later on proposed the Integer Projected Fixed Point 
(IPFP) algorithm [19], where they optimize the quadratic 
objective using a gradient-type algorithm interleaved with 
projection onto the discrete constraint set using e.g. the 
hungarian method. Lee et al. [15] tackle the graph match¬ 
ing problem using stochastic sampling, whereas Cho et 
al. [6] propose a reweighted random walk (RRWM) where 
the reweighting jumps are aiming at enforcing the match¬ 
ing constraints. Recently, Cho et al. [8] proposed a novel 
max pooling matching (MPM), where they tweak the power 
method to better cope with noise in the affinities. Although 
the algorithm comes without theoretical guarantees, it turns 
out to perform very well in practice, in particular, when one 
has a large number of outliers. Other second order methods 
include the work of Zhou and Torre [29] and Zaslavskiy 
et al. [25], where they propose deformable graph matching 
(DGM) and a path-following algorithm respectively. 

As second order methods are limited to pairwise similar¬ 
ities between two correspondences, higher order methods 
can integrate better geometric information. While they have 
higher computational complexity, they can cope better with 
geometric transformations such as scaling or other forms 
of noise. Compared to second order methods, higher or¬ 
der approaches are less well studied in literature due to the 
difficulty in handling the third order optimization problem. 

Duchenne et al. [10] formulated the hypergraph match¬ 
ing problem as a third order tensor optimization problem 
and proposed a higher order power method for solving the 
problem. This approach has shown significant improve¬ 
ments over second order ones which were proposed earlier. 
Zass and Shashua [26] introduced a probabilistic view to 
the problem, and proposed a Hypergraph Matching (HGM) 
algorithm. Their idea is to marginalize the tensor to a vec¬ 
tor and solve again a lower dimensional problem. Cher- 
tok and Keller [4] extended this idea and marginalized the 
tensor to a matrix, resulting in a second order matching 
problem which is solved by spectral methods. Both meth¬ 
ods are based on tensor marginalization, which leads to 


a loss of information. Moreover, they cannot effectively 
handle the one-to-one matching constraint during the iter¬ 
ations which is only considered at the final discretization 
step. Jungminlee et al. [16] extended the reweighted ran¬ 
dom walk approach of [15] to the third order setting. Their 
algorithm aims at enforcing the matching constraint via a 
hi-stochastic normalization scheme done for each iteration. 
In [27], Zeng et al. proposed to use pseudo-boolean opti¬ 
mization [1] for 3D surface matching, where higher order 
terms are decomposed into second order terms and then a 
quadratic pseudo-boolean optimization (QPBO) [13] algo¬ 
rithm is used to solve the problem. 

In this paper, we propose an algorithmic framework for 
solving the hypergraph matching problem based on a ten¬ 
sor block coordinate scheme. The key idea of our frame¬ 
work is to use a multilinear reformulation of the original 
objective function. In particular, we solve the third order 
graph matching problem by solving an equivalent fourth or¬ 
der one. Based on this reformulation, we derive two new hy¬ 
pergraph matching algorithms. We can guarantee for both 
algorithms monotonic ascent in the third order matching 
score. In the experiments we show that our new algorithms 
outperform previous work both in terms of achieving better 
matching score and accuracy. This holds in particular for 
very challenging settings where one has a high number of 
outliers and other forms of noise. 

2. Hypergraph Matching Formulation 

We formulate the correspondence problem as a hyper¬ 
graph matching problem. Hypergraphs allow modeling of 
relations which are not only pairwise as in graphs but in¬ 
volve groups of vertices. In this paper, we consider 3- 
uniform hypergraphs, that is each hyperedge contains 3 ver¬ 
tices. /c-uniform hypergraphs can be modelled as k-th order 
tensors which is the point of view we adopt in this paper. 

Given two attributed point sets G = {V, A) and G' = 
with Til = \V\ < n 2 = \V'\, the matching prob¬ 
lem can be formulated as finding a subset T' in the set of 
correspondences V x V' . The subset T’ can be represented 
by the binary assignment matrix X G where 

Xij = lif Vi e V matches Vj G V' and Xij = 0 else. A 
one-to-one matching X is an element of the set 

ni n2 

M = {X€ {0,1^""^ I < 1, ^Xij = 1}, 

i=l j=l 

that is we consider matchings which assign each vertex of V 
to exactly one vertex in V'. We further assume that we have 
a function : (V x V')^ M+ which assigns to each 
chosen triple C V and C V 

its similarity weight Finally, the score 


S of a. matching X can then be defined as 

ni n2 

E E , (*2 ,i2), (is Js) ^i2j2 ^isjs * 

ii,i2,is ji,j2,j3 

In order to facilitate notation we introduce a linear ordering 
inV X V' and thus can rewrite A as a vector x G {0, 
with n = nin 2 and becomes a tensor in so that 

the matching score, M, is finally written as 

n 

(1) 

i,j,k=l 

An m-th order tensor Q G - is called symmetric 
if its entries invariant under any permutation of 

their indices , im}- Note that is the sum of the 

componentwise product of the tensor X^ with the symmet¬ 
ric tensor Gfjj^ = x^x^x/^. Thus any non-symmetric part of 
X^ is “averaged” out and we can without loss of generality 
assume that X^ is a symmetric tensor^ In principle, one 
can integrate unary terms on the diagonal Xfa, i = 1,..., n 
and pairwise potentials Xfjj, i ^ j- However, the tensors 
we consider in Section 5 have just terms of order 3, that is 
Xfjk = f) if i = j, i = k or j = k. 

3. Mathematical Foundations for the Tensor- 
Block-Coordinate Ascent Method 

In this section we derive the basis for our tensor block 
coordinate ascent method to optimize directly over the 
set M of assignment matrices. The general idea is to op¬ 
timize instead of a homogeneous polynomial function the 
associated multilinear form. 

3.1. Lifting the Tensor and Multilinear Forms 

Instead of optimizing the third order problem in Equa¬ 
tion (1), we define a corresponding fourth order problem by 
lifting the third order tensor X^ to a fourth order tensor 

^ijki = ^ijk + ^iji + ^iki + ^jki' ( 2 ) 

This might be counter-intuitive at first in particular as pre¬ 
vious work [26, 4] has done the opposite way by reducing 
the third order problem to a second order problem. How¬ 
ever, lifting the tensor from a third order tensor to a fourth 
order tensor allows us to use structure of even order tensor 
which is not present for tensors of odd order. In particular, 
the score function is not convex. 

Lemma 3.1 Let be defined as in Equation (1). If is 

not constant zero, then 5'^ : ^ M is not convex. 

Tf is not already symmetric, then one can define E^ijk = 
f! E^ees K(i)a(3)a{ky where ©3 is the set of per- 

mutations of three elements. 




Proof: The Hessian HS^ of satisfies 

n 

HS^{x)jk = 6 ^ '^1 < 3, k< n, 

i=l 

for every x G i.e. iTS'^(x) = 6F^(x, •, •). Now, if 5'^ is 
convex, then 

0 < {y,HSH^)y) = 6F^^,y,y) Vx,y G R”. 

It follows that 

0 < F{x, y, y) and 0 < F(-x, y, y) = -F(x, y, y) 

for every x, y G In particular, for x = y we get 

0 = F(x, X, x) = S^{x.) for every x G □ 

The convexity is crucial to derive our block coordinate 
ascent scheme. As mentioned earlier, we do not work with 
the score function directly but with the multilinear form as¬ 
sociated to it. 

Definition 3.2 (Multilinear Form) The multilinear form 
X ... X ^ M associated to an m-th order 
tensor is given by: 

n 

F"(xi,...,x'") = y F™ixi --xr, 

il 

and the score function M associated to is 

defined by 5'"^(x) = F^(x,..., x). 

For simplicity, the superscript in the multilinear form 
can be omitted when there is no ambiguity. For ex¬ 
ample, F(x, y,z,t) can be used interchangeably with 
F^(x, y, z, t) in the same context since the number of vari¬ 
ables of this form is four which already implies a fourth or¬ 
der tensor. Also, we write F^(x, x, x, •) to denote the vector 
in such that 

n 

/^(x,x,x, ■)i = y F'^jkiXiXjXk, Vl<l<n. 

i,j,k=l 

Similarly, we write F^(x, x, •,•) to denote the matrix in 

I^nxn 

n 

F^{x,x,-,-)ki = y Vl<k,l<n. 

Note that if is associated to a symmetric tensor, then 
the position of the dots does not matter. For example, if 
is symmetric then one can check that /^(x, x, x, •) = 
F^(x, X, •,x) = F^(x, •,x, x) = F^(-,x, X, x). Similar 
properties hold for multilinear forms of other orders. 


It might seem at first that the lift from third to fourth 
order defined in (2) should lead to a huge computational 
overhead. However, 

n n 

S'^(x) = 4F^(x,x,x) y Xi = 4S'^(x) y Xi (3) 

and for all x G M, thus we have the fol¬ 

lowing equivalence of the optimization problems, 

max (x, x, x, x) = max F^ (x, x, x) (4) 

xGM xGM 

3.2. Convex Score Functions and Equivalence of 
Optimization Problems 

The use of multilinear forms corresponding to convex 
score functions is crucial for the proposed algorithms. The 
following lemma shows that even if we optimize the multi¬ 
linear form F^ instead of we get ascent in 5'"^. We can 
either fix all but one argument or all but two arguments, 
which leads to the two variants of the algorithm in Sec. 4. 

Lemma 3.3 Let F^ be a fourth order symmetric tensor. If 
5'^ : ^ M is convex, then for all x, y, z, t G 

1. /^(x,x,y,y) < max /^(u, u, u, u), 

uG{x,y} 

2. F^(x, y,z,t) < max F^(u, u, u, u). 

uG{x,y,z,t} 

Proof: First, we prove a characterization for the convexity 
of 5'"^. The gradient and the Hessian of 5'^ are given by 

V5"(x) = 4F^(x,x,x,.), 

HS‘^{x) = 12F^(x,x, 

where we used the symmetry of F^. 5'^ is convex if and 
only if iTS'^(x) is positive semi-definite for all x G 
which is equivalent to 

{y,HS\x)y) = \2F^{x,x,y,y) > 0,Vx,y G R”. (5) 

1. By (5) and the multilinearity of F, we have 

0 < F(x-y,x-y,x + y,x + y) (6) 

= F(x, X, X, x) + F(y, y, y, y) - 2F(x, x, y, y) 

for every x, y G It follows that 

2F(x,x,y,y) < F(x, x, x, x) + F(y, y, y, y) 

< 2 max F^(u,u,u,u). 

ue{x,y} 

for all X, y G 


n 


2. Similarly, for all x, y, z, t G we have 


and the proof is done. 


0 < F(x + y,x + y,z - t,z - t) 

= F(x, X, z, z) + F(x, X, t, t) + F(y, y, z, z) 

+2F(x,y,z,z) (7) 

+2F(x, y, t, t) - 2F(x, x, z, t) 

-2F(y, y, z, t) - 4F(x, y, z, t). 

Switching the variables from (x, y, z, t) to (z, t, x, y) 

and applying the same inequality, we get 

0 < F(x-y,x-y,z + t,z + t) 

= F(x, X, z, z) + F(x, X, t, t) + F(y, y, z, z) 
+^(y,y,t,t) - 2F(x,y,z,z) (8) 

-2F(x, y, t, t) + 2F(x, x, z, t) 

+2F(y, y, z, t) - 4F(x, y, z, t). 

Summing up inequalities (7) and (8) we obtain: 

4F(x,y,z,t) < F(x,x,z,z)+F(x,x,t,t) (9) 

+^(y,y,z,z) +F(y,y,t,t). 

Finally, applying the first result finishes the proof. 

n 

The following theorem shows that the optimization of the 
multilinear form and the score function are equiva¬ 
lent in the sense that there exists a globally optimal solution 
of the first problem which is also globally optimal for the 
second problem. 


As we cannot guarantee that 5'^ is convex for our chosen 
affinity tensor, we propose a modification making it convex. 
Importantly, this modification turns out to be constant on the 
set of matchings M. 

Proposition 3.5 Let be a fourth order symmetric ten¬ 
sor. Then for any a > where 11*^^112 • = 

frfffFTrtfy, the function 

S'^(x) := 5^(x) +Q; ||x ||2 
is convex on and, for any M, we have 
5'^(x) = 5'^(x) Fan\. 

Proof: The gradient and the Hessian of can be com¬ 
puted as: 

V5'^(x) = 4F^(x, X, X, •)+4a||x||2X, 

= 12F^(x, X, •, •) + 8axx^ + 4a ||x||2/, 

where / is the identity matrix. is convex if and only if 

(y,775^(x)y) >0 Vx,y 

This is equivalent to 

12/^(x, X, y, y) + Sa (x, y)^ + 4a ||x ||2 ||yII 2 > 0 (14) 


Theorem 3.4 Let be a fourth order symmetric tensor 
and suppose that M is convex. Then it holds for 

any compact constraint set D C BF, 

maxF^(x, X, X, x) = max /^(x,x,y,y) (10) 

= max F^(x, y,z,t) 

x,y,z,tG-D 

Proof: For any compact set D it holds: 

maxF^(x, X, X, x) < max F^(x, x, y,y) (11) 
xGD x,yGD 

< max F^(x, y,z,t). 

x,y,z,tG-D 

However, the second inequality in Lemma 3.3 shows 


for every x, y G B^. By the Cauchy-Schwarz inequality we 
have 


|F4(x,x,y,y)| = 




< 


Z] (Fjkif Z 






= ll-^lls 11^112 Ilyll 2 > 


for every x, y G B^. It follows that for a > 3 ||‘^^|| 2 ’ 
inequality (14) is true for any x, y G B^. Finally, the 
second statement of the proposition follows from the fact 
that ||x ||2 = rii for any x G M. □ 


/^(x,y,z,t) < max /^(u,u,u,u) (12) 

uG{x,y,z,t} 

which leads to 

max F^(x, y, z, t) < maxF^(x, X, X, x), (13) 

x,y,z,tG-D xG-D 


One of the key aspects of the algorithms we derive in 
the next section is that we optimize the multilinear form 
instead of the function 5'"^. This requires that we are able 
to extend the modification resp. ||x ||2 to a multilinear 
form. 







Proposition 3.6 The symmetric tensor Q'^ G 

with corresponding symmetric multilinear form defined as 


4. Tensor Block Coordinate Ascent for Hyper¬ 
graph Matching 


G'‘(x,y,z,t) = 


(x,y){z,t) + (x,z)(y,t) + (x,t){y,z) 


satisfies G'^{x,x,x,x) = ||x|| 2 . 

Proof: The proof requires two parts: 

1. There exists a fourth order symmetric tensor such 
that its multilinear form is given as above. 

2. The multilinear form associated to satisfies 

G'^(x,x,x,x) = ||x|| 2 . 

Let G be defined as follow: 


Sijki — ^ 


1 

II 

II 

II 

1/3 

if i = j fz k = 1 

1/3 

if i = k ^ j = 1 

1/3 

II 

II 

0 

else, 


VI < j, kj<n. 


The multilinear form associated to is then computed as 

n 

G'‘(x,y,z,t) = 

n 1 ^ 

= Ex«y«zA + - ^ XiyiZ.t, 

i=l ^ i,j = l 

^ n 1 ^ 

+ 3 H XjyjZjt^ + ^Y 

i,j=l i,j=^ 

= ^ ((x, y) (z, t) + (x, z) (y, t) + (x, t) (y, z)) 

As a result, we have G^(x, x, x, x) = ||x|| 2 . □ 

Thus in the algorithm we optimize finally 

/^(x, y, z, t) = /^(x, y, z, t) + a G^{x, y, z, t). 

Discussion: In [21] they proposed a general convexification 
strategy for arbitrary score functions where, similar to our 
modification, the added term is constant on the set of assign¬ 
ment matrices M. However, as the added term is inhomo¬ 
geneous, it cannot be extended to a symmetric multilinear 
form and thus does not work in our framework. In second 
order graph matching also several methods use convexified 
score functions in various ways [25, 28, 29]. However, for 
none of these methods it is obvious how to extend it to a 
third order approach for hypergraph matching. 


The key idea of both algorithms is that instead of directly 
optimizing the score function we optimize the corre¬ 
sponding multilinear form. Lemma 3.3 allows us then to 
connect the latter problem to the problem we are actually 
interested in. In both variants we directly optimize over the 
discrete set M of possible matchings, that is, there is no re¬ 
laxation involved. Moreover, we show mono tonic ascent for 
both methods. In most cases we achieve higher scores than 
all other higher order methods, which is reflected in signifi¬ 
cantly better performance in the experiments in Section 5. 

The original graph matching problem is to maximize the 
score function over all assignment matrices in M. This 
combinatorial optimization problem is NP-hard. As stated 
above, instead of working with the original score function 
we use the lifted tensor and thus the lifted score func¬ 
tion S^. As we cannot guarantee that is convex, we 
optimize finally resp. the associated multilinear form 

. However, we emphasize that even though we work with 
the lifted objects, both algorithms prove monotonic ascent 
with respect to the original score function over the set 
M. The central element of both algorithms is Lemma 3.3 
and the idea to optimize the multilinear form, instead of the 
score function directly, combined with the fact that for as¬ 
signment matrices both modifications (lifting and convexifi¬ 
cation) are constant and thus do not change the problem. In 
order to simplify the notation in the algorithms we discard 
the superscript and write instead of as the algorithms 
only involve fourth order tensors. 

4.1. Tensor Block Coordinate Ascent via Linear As¬ 
signment Problems 

The first algorithm uses a block coordinate ascent 
scheme to optimize the multilinear form F ^. In particular, 
we solve the following optimization problem 

max F^(x,y,z,t). 

x,y,z,tGM 

Fixing all but one argument and maximizing that over all 
assignment matrices leads to a linear assignment problem 
which can be solved globally optimal by the Hungarian al¬ 
gorithm [14, 2] - this is used for steps l)-4) in Algorithm 1. 
However, optimization of the multilinear form is not neces¬ 
sarily equivalent to optimizing the score function. This is 
the reason why we use Lemma 3.3 in step 5) to come back 
to the original problem. The following theorem summarizes 
the properties of Algorithm 1 . 

Theorem 4.1 Let be a fourth order symmetric tensor. 
Then the following holds for Algorithm 7.* 

7. The sequence 7^q,(x^, y^, z^, t^) for k = 1, 2,... A 
strictly monotonically increasing or terminates. 




and 


Algorithm 1: BCAGM 
Input: Lifted affinity tensor a = 3 || 

(x^,y^,z^, t^) G M X M X M X M, /c = 0, m = 0 

Output: X* G M 
Repeat 


1. 

^fc+i ^ 

argmaXxeM^’aXy'' 


2. 

yfe+l ^ 

: argmaXygMf’a(x''+^ 

,y,z\C) 

3. 

^k+1 ^ 

arg maXz£Mf’a(x'‘+\ 

++i,z,t''; 

4. 

II 

+ 

arg maxteMf’a(x''+\ 

■{V/v+l ^/c+1 

y 5 

5. 

ifF„(x 

/c+1 y/C+1 ~/C+l .j-/C+l^ 

II 


then 




- = argmax Fq,(u, u, u, u) 

- ifF„(u™+\u'"+\u™+\u'"+i) = 

Fq,(x^+^, y^+^, z^+^, then return 

_~ 

- m = m -\-1 

else x^^^ = x^^^, = y^^^, z^^^ = z^^^, 


end 

6. k = k -\-1 


2. The sequence of scores 5'^(u"^) for m = 1, 2,... A 
strictly monotonically increasing or terminates. For 
every m, E M is a valid assignment matrix. 

3. The sequence of original third order scores S^{u'^) 
for m = 1,2,... is strictly monotonically increasing 
or terminates. 

4. The algorithm terminates after a finite number of iter¬ 
ations. 

Proof: From the definition of steps 1) — 4) in Algorithm 1, 
we get 

Fa,k ■■= Fa{x'",y'",'z'",C) 

< Fa{±’^+\y^z^V) 

< F„(x'=+^y'=+^z^t'') 

< F„(x'=+\y*+\z*+\t'') 

< F„(x'=+i,y'=+i,z'=+i,P+i) =; F„,fe+i. 

Either F^^k+i > ^a,/c in which case 

+ l /c + 1 ~/c+l /c + 1 ~/c + l ^/c + 1 ^/c+1 

X —X ,y —y ,z —z ,t — t 


F„(x'=+i,y'=+\z'=+\t''+i) > F„(x^y^z^t'=), 

or Fa,k+i = Fa,k and the algorithm enters step 5. Since, 
by Proposition 3.5, 5'^ is convex for the chosen value of a, 
applying Lemma 3.3 we get 

Fa,k+i < max _ F„(v,v,v,v) 

v^l^x^+i 1 

= F„(u'"+\u™+\u'"+\u™+i) 

= 54(u™+1). 

If the inequality is an equality, then the termination condi¬ 
tion of the algorithm is met. Otherwise, we get 

Fa,k+1 < F„(u™+I,u™+I,u'"+I,u™+1) 

= 54(u™+ 1) =F„(x'=+\y'=+\z'=+\V+i). 

This proves the first statement of the theorem. 

From 

= F.(x^+\y^+\z^+\t^+i) 

it follows that 5'^(u’^), m = 1, 2,... is a subsequence of 
Fq,(x^, y^, z^, t^),/c = 1,2,... and thus it holds either 
in which case the algorithm ter¬ 
minates or However, by Equation 

(3.5), the additional term which has been added to to get 
a convex function is constant on M, that is 

5'^(x) = S^{x.) + anl Vx G M. 

It follows that either and the algo¬ 
rithm terminates or S^{u^) < which proves the 

second part of the theorem. 

By Equation (3), we have 5'^(x) = 4ni 5'^(x) for any 
X G M. Thus the statements made for directly translate 
into corresponding statements for the original third order 
score S^. This proves the penultimate statement. 

Einally, the algorithm yields a strictly monotonically in¬ 
creasing sequence 5'^(u"^), m = 1, 2,... or it terminates. 
Since there is only a finite number of possible assignment 
matrices, the sequence has to terminate after a finite number 
of steps. n 

We would like to note that all statements of Theorem 4.1 
remain valid for a = 0 if 5'^ is convex. This is also the 
motivation why we run the algorithm first with a = 0 until 
we get no further ascent and only then we set a = 3|| J^|| 2 . 
It turns out that in the experiments often the phase of the al¬ 
gorithm with a = 0 leads automatically to a homogeneous 
solution and no further improvement is achieved when set¬ 
ting a = 3IIII 2 . This suggests that the constructed score 
functions in the experiments are already convex or at least 
close to being convex. 





4.2. Tensor Block Coordinate Ascent via a Sequence 
of Quadratic Assignment Problems 

The second algorithm uses a block coordinate ascent 
scheme to optimize the multilinear form Fcx where now two 
arguments are fixed and one maximizes over the other two. 
The resulting problem is equivalent to a quadratic assign¬ 
ment problem (QAP) which is known to be NP-hard. Thus a 
globally optimal solution as for the linear assignment prob¬ 
lem in Algorithm 1 is out of reach. Instead we require a 
sub-algorithm which delivers monotonic ascent for the 
QAP 

max (x, Ax), (15) 

xGM 

where A e is a nonnegative symmetric matrix. As 

in Algorithm 1 , we go back to the optimization of the score 
function in step 3) using Lemma 3.3. The following theo¬ 
rem summarizes the properties of Algorithm 2. 


Algorithm 2: BCAGM-T^ 

Input: Lifted affinity tensor a = 3 || 

(x^,y^) G M X M, k = 0, m = 0, 
z = x^) is an algorithm for the QAP in (15) 

which provides monotonic ascent, that is 
(z, Az) > (x^, Ax^) 

Output: X* G M 
Repeat 

1. x'=+i = «'(F„{.,.,yGy'=), 

2. y''+i = «'(F„(5c*+i,x''+i,-,-),y'') 

3. ifF„(5c'=+i,x'=+i,y'=+i,y''+i) = F„(x^x^y^y'') 

then 

_ = argmax Fq,(u,u, u, u) 

uG{x^+^ ,y^+^ } 

- if F„(x''+\x''+\y*+\y*^+i) = 

Fa(u"^+\ u™+i, u^+i, u™+i) then return 

- x^~^^ = 

- m = m -h 1 

else x^~^^ = x^~^^ ~ 

end 

4. k = k 1 


Theorem 4.2 Let be a fourth order symmetric tensor 
and let be an algorithm for the QAP which yields mono¬ 
tonic ascent. Then the following holds for Algorithm 2: 

1. The sequence Fq,(x^, x^, y^, y^) for k = 1, 2,... A 
strictly monotonically increasing or terminates. 


2. The sequence of scores 5'^(u"^) for m = 1, 2,... A 
strictly monotonically increasing or terminates. For 
every m, ^ M is a valid assignment matrix. 

3. The sequence of original third order scores S^{u'^) 
for m = 1,2,... is strictly monotonically increasing 
or terminates. 

4. The algorithm terminates after a finite number of iter¬ 
ations. 


Proof: From the definition of steps 1) — 2) in Algorithm 2, 
we get 

:= F„(x^x^y^y'=) 

< F„(x''+^x'=+^y^y'') 

< F„(x''+\x'=+\y''+\y''+i) =: F„,fc+i. 

Either Fo,^k+i > ^a,/c in which case 

X^+l = x^+l y^ + 1 = y^ + 1 


and 


F„ (x'^+i, x'^+i, y'^+i, y'^+i) > F„ (x^ x^ y^ y'^), 

OX Fo,^k+i = and the algorithm enters step 3. Since, 
by Proposition 3.5, S% is convex for the chosen value of a, 
applying Lemma 3.3 we get 


-^ctj/c + l ^ 


max Fq,(v,v,v,v) 

vG{x^+^ 


= F„(u'"+\u™+\u'"+\u™+i) 


If the inequality is an equality, then the termination condi¬ 
tion of the algorithm is met. Otherwise, we get 

K,k+i < F„(u™+\u”*+\u™+\u”*+i) 

= =F„(x'=+i,x''+i,y''+i,y'=+i). 


This proves the first statement of the theorem. 
From 


5'4(u”*+ 1) = F„(u”*+\u™+\u”*+\u™+i) 

= F„(x'=+\x''+i,y'=+i,y'=+i) 

it follows that = 1, 2,... is a subsequence of 

Fq,(x^, x^, y^, y^), k = 1,2,... and thus it holds either 
5'^(u"^) = 5'^(u"^+^) in which case the algorithm termi¬ 
nates or S'^(u"^) < 5'^(u’^+^). However, by Equation 
(3.5), the additional term which has been added to to 
get a convex function is constant on M, that is 

5'^(x) = 5'^(x) + anf Vx G M. 





It follows that either and the algo¬ 
rithm terminates or >5^ (u"^) < 5'^ which proves the 

second part of the theorem. 

By Equation (3), we have 5'^(x) = 4ni 5'^(x) for any 
X G M. Thus the statements made for 5'^ directly translate 
into corresponding statements for the original third order 
score S^. This proves the penultimate statement. 

Finally, the algorithm yields a strictly monotonically 
increasing sequence m = 1, 2,... or it terminates. 

Since there is only a finite number of possible assignment 
matrices, the sequence has to terminate after a finite number 
of steps. n 

In analogy to Theorem 4.1 all statements of Theorem 4.2 
remain valid for a = 0 if is convex. Thus we use the 
same initialization strategy with a = 0 as described above. 
There are several methods available which we could use for 
the sub-routine in the algorithm. We decided to use the 
recent max pooling algorithm [8] and the IPFP algorithm 
[19] and use the Hungarian algorithm to turn their output 
into a valid assignment matrix in M. It turns out that the 
combination of our tensor block coordinate ascent scheme 
using their algorithms as a sub-routine yields very good per¬ 
formance on all datasets. 

5. Experiments 

We evaluate our hypergraph matching algorithms on 
standard synthetic benchmarks and realistic image datasets 
by comparing them with state-of-the-art third order and sec¬ 
ond order approaches. In particular, we use the following 
third order approaches: Tensor Matching (TM) [10], Hyper¬ 
graph Matching via Reweighted Random Walks (RRWHM) 
[16], Hypergraph Matching (HGM) [26]. For second 
order approaches. Max Pooling Matching (MPM) [8], 
Reweighted Random Walks for Graph Matching (RRWM) 
[6], Integer Projected Fixed Point (IPFP) [19], and Spec¬ 
tral Matching (SM) [17] are used. We denote our tensor 
block coordinate ascent methods as BCAGM from Algo¬ 
rithm 1 which uses the Hungarian algorithm as subroutine, 
and BCAGM-fMP and BCAGM-fIPFP from Algorithm 2 
which uses MPM [8] and IPFP [19] respectively as subrou¬ 
tine. MPM has recently outperformed other second order 
methods in the presence of a large number of outliers with¬ 
out deformation. Thus, it serves as a very competitive can¬ 
didate with third order approaches in this setting. For all the 
algorithms, we use the authors’ original implementation. 

In all experiments below, we use the all-ones vector as 
starting point for all the methods. Moreover, we apply the 
Hungarian algorithm to the output of every algorithm to turn 
it into a proper matching. 

Generation of affinity tensor/matrix: To build the 
affinity tensor for the third order algorithms, we follow the 


approach of Duchenne et al. [10]: for each triple of points 
in the same image we compute a feature vector / based 
on the angles of the triangle formed by those three points. 
Let fji,j 2 ,j 3 denote two feature vectors for two 

triples {ii,i 2 ,h) and (ji,^ 2 , 73 ) in two corresponding im¬ 
ages. Then we compute the third order affinity tensor as: 

^(ildl),(i2,j2),ii3,j3) ~ exp ( — 7 II/ii,22 43 ~ /jl ,72 JS 11 2 ) 

(16) 

where 7 is the inverse of the mean of all squared distances. 
As shown by Duchenne et al. [10] these affinities increase 
the geometric invariance compared to second order models 
making it more robust to transformations, like translations, 
rotations and scalings. This partially explains why third or¬ 
der approaches are preferred to second order ones in this 
regard. As the computation of the full third order affinity 
tensor is infeasible in practice, we adopt the sampling strat¬ 
egy proposed by [10, 16]. 

The affinity matrix for all the second order methods are 
built by following [8, 6], where the Euclidean distance of 
two point pairs is used to compute the pairwise similarity: 

- 4i2lV^a) (17) 

where the Euclidean distance between two points 

A, ^2 ^ P, is distance between the points 71,72 ^ 
Q and as is a normalization parameter specified below. 

5.1. Synthetic Dataset 

This section presents a comparative evaluation on match¬ 
ing point sets P and Q in . We follow the approach in 
[8] for the creation of the point sets. That is, for P we sam¬ 
ple Tiin inlier points from a Gaussian distribution A/’(0,1). 
The points in Q are generated by adding Gaussian deforma¬ 
tion noise A/’(0, cr^) to the point set P. Furthermore, Uout 
outlier points are added to Q sampled from A/’(0,1). We 
test all matching algorithms under different changes in the 
data: outliers, deformation and scaling (i.e. we multiply the 
coordinates of the points in Q by a constant factor). 

In the outlier setting, we perform two test cases as fol¬ 
lows. In the first test case, we vary the number of outliers 
from 0 to very challenging 200, which is 20-times more than 
the number of inliers riin, while cr is set to 0.01 and no scal¬ 
ing is used. In the second case, we set a to 0.03 and scale all 
the points in Q by a factor of 1.5, making the problem much 
more difficult. This test shows that third order methods are 
very robust against scaling compared to second order meth¬ 
ods. In all the plots in this section, each quantitative re¬ 
sult was obtained by averaging over 100 trials. The accu¬ 
racy is computed as the ratio between the number of correct 
matches and the number of inlier points. We use as = 0.5 
in the affinity construction for the second order methods as 
suggested by [8]. The results in Figure 1 show that our al¬ 
gorithms are robust w.r.t. outliers even if deformation and 
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Figure 2: Robustness of matching algorithms w.r.t. defor¬ 
mation. (riin = 20, flout = 0, scale = 1.) 


scaling are involved. When the deformation is negligible 
as shown in Figure la, MPM and BCAGM-fMP perform 
best followed by BCAGM, BCAM-fIPFP and RRWHM. In 
a more realistic setting as shown in Figure lb, where out¬ 
lier, deformation and scaling are present in the data, our 
algorithms outperform all other methods. One can state that 
BCAGM-fMP transfers the robustness of MPM to a third 
order method and thus is able to deal with scaling which is 
difficult for second order methods. 

In the deformation setting, the deformation noise a is 
varied from 0 to challenging 0.4 while the number of inliers 
riin is fixed to 20. There are no outlier and scale changes 
in this setting. This type of experiment has been used in 
previous works [8, 16]. where a has been varied from 0 to 
0.2. The result in Figure 2 shows that our algorithms are 
quite competitive with RRWHM which has been shown to 
be robust to deformation [16]. It is interesting to note that 
while MPM is very robust against outliers (without scal¬ 
ing) it is outperformed if the deformation is significantly in¬ 
creased (without having outliers). However, note that even 
though our BCAGM-i-MP uses MPM as a subroutine, it is 
not affected by this slight weakness. Considering all ex¬ 
periments our BCAGM-i-MP and to a slightly weaker extent 
BCAGM and BCAGM-fIPFP outperform all other methods 
in terms of robustness to outliers, deformation and scaling. 
Runtime: Figure 1 shows that the running time of our al¬ 
gorithms is competitive in comparison to other higher order 
and even second order methods. In particular, BCAGM is 
always among the methods with lowest running time (more 
results can be found in the appendix). 

5.2. CMU House Dataset 

The CMU house dataset has been widely used in previ¬ 
ous work [6, 16,10, 29] to evaluate matching algorithms. In 
this dataset, 30 landmark points are manually tracked over a 
sequence of 111 images, which are taken from the same ob¬ 
ject under different view points. In this experiment, “base¬ 
line” denotes the distance of the frames in the sequence and 
thus correlates well with the difficulty to match the corre¬ 
sponding frames. 

We matched all possible image pairs with “baseline” 
of 10, 20,30,..., 100 frames and computed the average 



(a) 34 pts vs 44 pts, 10 outliers 


(b) TM 10/34(1715.0) 





(c) HGM 9/34 (614.7) 


(d) RRWHM 28/34 (5230.5) 



(e) BCAGM 28/34 (5298.8) 


(f) BCAGM+MP 34/34 (5377.3) 


Figure 4: Car dataset: the number of correct matches and 
the objective score are reported. (Best viewed in color.) 





(e) BCAGM 19/23 (4016.2) 


(f) BCAGM+MP 16/23 (4133.6) 


Figure 5: Motorbike dataset: the number of correct matches 
and the objective score are reported. (Best viewed in color.) 


matching accuracy for each algorithm. The algorithms are 
evaluated in three settings. In the first experiment, we match 
30 points to 30 points. Then we make the problem signifi¬ 
cantly harder by randomly removing points from one image 
motivated by a scenario where one has background clutter 
in an image and thus not all points can be matched. This 
results in two matching experiments, namely 10 points to 
30 points, and 20 points to 30 points. For the choice of ds 
in the affinity tensor for second order methods, we follow 
[6, 29] by setting ds = 2500. 

The experimental results are shown in Figure 3. While 
most algorithms perform rather well on the 30 to 30 task, 
our methods perform significantly better than all other 
methods in the more difficult tasks, thus showing as for the 
synthetic datasets that our methods are quite robust to dif¬ 
ferent kind of noise in the matching problem. 
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Figure 1: Matching point sets in each row shows the accuracy, matching score and running time of all algorithms. The 
number of outliers has been varied from 0 to 200. (a) Increasing number of outliers with slight deformation and no scaling, 
(b) Increasing number of outliers with larger deformation and scaling. (Best viewed in color.) 




(a) Car (b) Motorbike 

Figure 6: Evaluation of higher order GM algorithms on the 
car and motorbike dataset. (Best viewed in color.) 

5.3. Car Motorbike Dataset 

In this experiment, we compare our algorithms with 
other third order approaches on the car and motorbike 
dataset introduced in [18]. The dataset contains 30 resp. 
20 pairs of car resp. motorbike images. Each image con¬ 
tains a number of inliers and outliers from the cluttered 
background. Ground-truth correspondences are given for 


the inlier points in both images. Figures 4 and 5 show 
some examples of matching results. To test the algorithms 
against noise, we kept all the inlier points in both images 
and randomly added a number of outliers from the cluttered 
background to the second image (between 0 to 15). For 
each number of outliers, we test all algorithms on all the 
image pairs of each dataset and report the average perfor¬ 
mance in Figure 6. It shows that our methods achieve bet¬ 
ter overall performance than other third order approaches 
in both datasets. However, we expect that the performance 
can be significantly improved if also other features are in¬ 
tegrated rather than just relying on geometric information 
[28, 29], or graph learning methods can be also be inte¬ 
grated [5, 3, 18, 22, 24]. 

6. Conclusion and Outlook 

We have presented a new optimization framework for 
higher order graph matching. Compared to previous meth¬ 
ods for higher order matching we can guarantee monotonic 
ascent for our algorithms in the matching score on the set of 
assignment matrices. Our new algorithms achieve superior 
performance in terms of objective but also yield competi¬ 
tive or significantly better matching accuracy. This is par¬ 
ticularly true for large number of outliers and other forms 
of noise. Moreover, both algorithms are also competitive in 




































(a) An input pair: 10 pts vs 30 pts, baseline = 50 


(b) MPM 4/10 (15.6) 


(c) TM 5/10(18.4) 



(d) RRWHM 7/10 (26.4) (e) BCAGM 10/10 (43.1) (f) BCAGM+MP 10/10 (43.1) 
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Figure 3: CMU house dataset: the first two rows show the matching output of several algorithms on an example pair of 
images where the baseline is 50. The yellow/red lines indicate correct/incorrect matches. The third row shows the average 
performance of algorithms with different number of points in the first image. (Best viewed in color.) 


terms of runtime compared to other methods. 

An interesting line of future research is the use of glob¬ 
ally optimal solutions of relaxations of the hypergraph 
matching problem. We computed via the method of [12] 
the maximal singular vectors of the fourth order tensor and 
used these as initialization of Algorithm 1 . This lead to an 
improvement in matching score of up to 2% and in accuracy 
of up to 3%. We will explore this further in future work. 
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Appendix 

In this section, we show some additional experimental 
results to the main paper. 

A. Synthetic Dataset 

Figure A shows additional results in the outlier setting 
where the number of inliers was fixed to 10 while a and 
scale were set to different values. It is interesting to see 
that when there is no deformation and scaling, our algo¬ 
rithms together with RRWHM [16] and MPM [8] achieve 
an almost perfect result. However, our algorithms outper¬ 
form all other higher order approaches when deformation 
and scaling are slightly present. Compared to second order 
methods, our algorithms can take advantage of higher or¬ 
der features, therefore, achieve superior performance when 
transformations such as scaling are present. 

Figure B shows further results in the deformation setting, 
where the number of inliers was set to 30 and 40 accord¬ 
ingly, and no other form of noise was used. As we show 
the runtime for all the experiments, the result for 20 inliers 
points is repeated from the paper as well. One can observe 
from Figure B that our algorithms always stay competitive 
with other state-of-the-art higher order methods, in particu¬ 
lar RRWHM [16], even when the deformation is significant. 

B. CMU House Dataset 

Similar to the experiments done in Section 5.2, we eval¬ 
uate GM algorithms on two tasks where we match 15 points 
to 30 points and 25 points to 30 points in two correspond¬ 
ing images. For each task, we match all the possible im¬ 
age pairs and compute the average result for each baseline. 
The results in Figure C show that our algorithms achieve 
competitive or better results than other methods for all the 
baselines. 

C. Car Motorbike Dataset 

Figure D shows the running time of higher order meth¬ 
ods for the experiments done in Section 5.3 on the Car and 
Motorbike dataset. 
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Figure A: Matching point sets in (outliers test): The top row shows average accuracy while the middle row shows the 
average matching score and the bottom row shows average running time. The number of outliers was varied from 0 to 200 
with the interval 10. (a) Increasing number of outliers without deformation and scaling, (b) Increasing number of outliers 
with slight deformation and small scaling, (c) Increasing number of outliers with slight deformation and large scaling. (Best 
viewed in color.) 
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Figure B: Matching point sets in (deformation test): The top row shows average accuracy while the middle row shows 
average objective score and the bottom row shows average running time, (a) 20 inliers (b) 30 inliers (c) 40 inliers (Best 
viewed in color.) 
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Figure C: CMU house dataset: The top row shows average matching accuracy while the middle row shows average objective 
score and the bottom row shows average running time. (Best viewed in color.) 
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Figure D: Car and Motorbike dataset: Running time of the third order methods. (Best viewed in color.) 


































